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INTRODUCTION 



The field of three-dimensional (3D) underwater acoustics has been 
developed in an attempt to describe the properties and characteristics of 
sound propagation in the heterogeneous ocean environment (Lynch and 
Chiu, 1989). At present, there are three major approaches used in studying 
the 3D underwater sound field. One of these, ray theory, is premised upon a 
high-frequency approximation method. This method gives a geometric- 
optics solution which involves ray-tracing in a spatially varying sound speed 
field. The Hamiltonian Acoustic Ray Tracing Program for the Ocean 
(HARPO) is a versatile numerical code developed by Jones et al. (1986) for the 
computation of 3D rays. The primary deficiency of this theory is that it cannot 
adequately address the behavior of low-frequency sound due to the neglect of 
sound dispersion and diffraction. 

The second approach uses a parabolic approximation to the acoustic wave 
equation, which was introduced by Tappert (1977). A 3D numerical parabolic 
equation (PE) model was developed by Lee et al. (1988) using an implicit finite 
difference scheme. Another 3D PE model was arrived at by Baer (1981) 
utilizing a split-step Fourier algorithm. In a recent study, analytic solutions to 
the 3D parabolic equation were obtained by Seigmann et al., (1990). These 
solutions are valuable for testing the accuracy of 3D numerical models. 

A third approach is the normal mode method. Pierce presented a three- 
dimensional version of this method in 1965. The normal mode method is 
based upon a separable solution to the wave equation. Pierce assumed an 
adiabatic acoustic environment which leads to the neglect of coupling 



1 



presented by Chiu and Ehret (1990). This later 3D normal mode model 
accounts for both horizontal refraction and model coupling. 

Any improvement in acoustic modeling must be quantified. We shall 
endeavor to test the accuracy of the 3-D coupled normal mode model of Chiu 
and Ehret (1990) in regard to horizontal refraction. 

To examine the accuracy of the Chiu-Ehret model, we compare the results 
of their model with those from analytic solutions to the parabolic equation 
arrived at by Seigmann, et al. (1990). The results of the Chiu-Ehret model and 
the analytic solutions of the parabolic equation are compared for two cases of 
horizontal variation in sound speed. As a means of comparing the results, we 
examine the slow variations of the complex pressure envelope function and 
transmission loss. 

Our simulated acoustic field has a range of 100 km with azimuth from 30° 
to 180°. Two variations in the horizontal sound speed field were used. In 
Case I the sound speed varies only with the azimuth angle while Case II varies 
both with the azimuth and radial. 

The normal mode theory used by Chiu and Ehret is described in detail in 
the Appendix. The appendix also includes a description of the PE 
approximation and the analytic solutions developed by Seigmann et al. 
(1990). 
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II. ANALYSIS AND METHOD OF MODEL COMPARISON 



The three-dimensional coupled normal mode model was applied to two 
sound speed fields for which exact analytic solutions to the PE approximation 
are available. These analytic PE solutions for vertically invariant sound speed 
were developed by Seigmann et al. (1990). Their development was based on 
expressing the "time-independent" acoustic pressure component as 



0pE{r,e,z) = Up{r,e)smi^YjZ^exp 
where 
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Cq = reference sound speed 
/ = acoustic frequency 
H = ocean depth. 

The reduced complex envelope function 
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where 
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«(r) = ^ {r = {r,0,z)) 

C(r) = sound speed as a function of position r. 



Substitution and separation of real and imaginary parts gives two coupled 
equations governing and 0: 



dAp Ap d^0 

dr kQr^ dd do 2kQr^ dO^ 
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(2.3) 



If one chooses 0{r,9) = a{r)6, equation (2.2) may be solved for the amplitude, 
i.e.. 





(2.4) 



Choosing a(r) and the functional dependence F results in a reduced 
envelope function (and ultimately acoustic pressure) while equation (2.3) can 
be solved for the index of refraction n{r) and thus the sound speed field C(r). 
Equation (2.3) is for a wide angle PE approximation. This wide angle form 
will be used throughout. To judge the accuracy of the Chiu-Ehret model in 
calculating horizontal refraction we compare the numerical coupled mode 
results with the analytic PE results for two sound speed fields or two 
functional definitions of (x{r). 
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A. INDEX OF REFRACTION 



1. Case I 

For Case 1 we choose a(r) = ocQkQ r, or 

0 = rG (2.5) 



where «q is an arbitrary constant. From equation (2.4) the amplitude is 

A^(r,G) = F(0-aglnr). 



Let the function F be multiplied by amplitude has the final form 



^/r,G) = /3gA'o(0-aolnr) 



(2.6) 



where j3 g is an arbitrary constant with units of length. The square of the index 
of refraction from (2.5) and (2.6) is 
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The resultant sound speed is shown in Figure 1 for ag = 0.005 and 
Cg = 1500 m/s. 

2. Case II 

For Case II we choose 

0 = a^k^e (2.8) 



and hence 

A^(r,0) = j3g/:g(0-ag^-gr). (2.9) 

The square of the index of refraction is now 
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The resulting sound speed is shown in Figure 2 for «q = 4x10 and Cq = 1500 
m/s. 

B. ENVELOPE FUNCTION 

We now have two analytic sound speed fields and the associated acoustic 
pressure fields derived from the parabolic equation approximation. 
Considering only j = 1, Seigmann et al., (1990) define pressure as 

1 



0PE{r, d,z) = 



ynkQVj 



2 ^2 '1^0''“] 

Mpg'® sin(7|z)exp 



ikQV 
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4 






( 2 . 11 ) 



For Case I, using the amplitude and phase expressed in (2.6) and (2.5), 
respectively, we obtain 

0pE{r,e,z) = (-^Ye^''^''^\okQ{e-ao\nr)e^^^^^^ 
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xexp 
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( 2 . 12 ) 



For Case II, using (2.8) and (2.9), we get 
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Figure 1. Case I Sound Speed Field for Uq = 0.005 and = 1500 m/s 
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(2.13) 
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0PE{r,e,z) = \^—Ye'^^ '^^Poko{e-aokor)e^^^o’' ^sin(yiz) 
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In terms of normal modes, pressure can be expressed as 

<*’NM('-,e,2)= (2-M) 

th 

where is the n normal mode (eigenmode), k^ is the corresponding 
horizontal wavenumber (eigenvalue) and li„ is the corresponding slowly 
varying modulation envelope. See the appendix for details. 

Since we are considering depth-invariant sound speed fields, an analytic 
solution is available for the eigenvalues and eigenmodes. The first mode is 




and thus, keeping only the contribution of the first mode, (2.14) becomes 

^NM 2) = (^' sin( 7iz). (2. 15) 

The normal mode model entails a numerical calculation of the slowly 
varying envelope function U^. Therefore, the quantification of model 
accuracy can be achieved by comparing the numerical l/„ to those derived 
from the analytic PE solutions. The analytic envelope function U.[(r,0) for 
Case I is derived by equating (2.12) to (2.15). The resulting expression is 
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With the index of refraction n given by (2.7), the range integral of k-^ can be 
expressed as 
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It follows that the analytic envelope function U-^(r,6) can be recast as 
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For Case II, the wave number integral, using (2.10), becomes 
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The analytic envelope function is then 
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For Case I we use the analytic envelope function (2.17) to define 
the initial condition at r = Tq for the numerical model run. The condition is 
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For Case II we use the analytic envelope function from (2.18) at r = Tq. The 
Case II initial condition is 



Ui{ro,9) 




'Poh{9 - ocQkQVQ) 




( 2 . 20 ) 
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0 = -7«o^o'o‘ 

o 
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III. MODEL RUNS 



Our computational domain is in a cylindrical coordinate system. The 
horizontal range is from 1 km to 100 km. The azimuth range is 30° to 150°. 
The depth is 4000 m. The domain is the same for both cases. 

We use a frequency of 50 Hz and source depth of 1000 m. The reference 
sound speed is 1500 m/s. We consider mode 1 only. The vertical boundary 
conditions for the problem are pressure release surface and rigid bottom. Two 
range and azimuth dependent sound speed fields as presented in the last 
chapter are used. For Case 1 sound speed ranges from 1481 m/s to 1496 m/s, 
varying only in azimuth with a gradient of 0.128 to 0.130 m/s/degree (see 
Table 1). For Case II sound speed varies in both range and azimuth (see Table 
1 again). Here sound speed goes from 1437 m/s to 1500 m/s with a radial 
gradient of 0.13 to 0.657 m/s/km and an azimuth gradient of 0.004 to 0.431 
m/s/ degree. Both selections of sound speed fields contain realistic gradients 
as observed in the ocean. 



TABLE 1. THE AZIMUTH AND RADIAL SOUND SPEED VARIATION 



Case 


Range 

r[km] 


Sound Speed [m/s] 


Azimuth Variation 

i-T'-sl 


Radial Variation 
dC m 






0=30° 


0=90° 


0=150° 


0=30° 


0=90° 


0=150° 


0=30° 


0=90° 


0=150° 


I 


1 


1496 


1488 


1481 


0.1303 


0.1289 


0.1275 


0 


0 


0 


50 


1496 


1488 


1481 


0.1303 


0.1289 


0.1275 


0 


0 


0 


100 


1496 


1488 


1481 


0.1303 


0.1289 


0.1275 


0 


0 


0 


II 


1 


1500 


1500 


1499 


0,0044 


0.0043 


0.0043 


0.1316 


0.3944 


0.6482 


50 


1493 


1480 


1468 


0.2217 


0.2178 


0.2142 


0.1310 


0.3848 


0.6221 


100 


1487 


1462 


1437 


0.4265 


0.4127 


0.4001 


0.1301 


0.3763 


0.5992 
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Eigenmodes and eigenvalues are calculated using finite difference 
approximations to (A. 8) in the appendix with a 40 m vertical grid spacing and 
a matrix eigenvalue solver. 

To obtain the numerical envelope function a first order differential 
equation ((A. 24) in the appendix) is integrated. The radial integration step 
size is 0.5 km. The 61 integration paths were separated by two degrees. 
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IV. RESULTS 



A. CASE I 

In comparing the numerical normal mode (NM) and analytic parabolic 
equation (PE) results we focus on the slowly varying modulating pressure 
envelope. The removal of the rapid oscillations, which have wave lengths of 
order 27 r/l<Q, makes the displays of results easier to interpret. 

Figures 3 and 4 show the amplitude and phase, respectively, of the 
envelope function calculated by the normal mode model using the Case I 
sound speed field. Figures 5 and 6 are the amplitude and phase, respectively, 
of the envelope based on the analytic solution to the parabolic equation. 
Figure 7 shows the relative difference, defined as 



The percent difference is everywhere less than 5%. The difference in NM and 
PE phase is shown in Figure 8. 

Transmission loss at a depth of 1000 m calculated by the NM model, 
TLj^^j, is shown in Figure 9. The magnitude of the analytic PE pressure, from 
(2.12), is 



Upg( analytic solution) - (numerical solution) 



Upp (analytic solution) 




(3.1) 



Substituting (2.6) in (3.1), we get 
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Figure 3. Amplitude of for Case I with a = 0.005 and P = 0.1 




Figure 4. Phase of for Case I with a = 0.005 and P = 0.1 
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Figure 5. Amplitude of Upg for Case I with a = 0.005 and = 0.1 




Figure 6. Phase of Upg for Case I with a = 0.005 and /3 = 0.1 
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Figure 7. Amplitude Difference in Percent for Case I with a = 0.005 and 

P = 0.1 




Figure 8. Phase Difference in Degree for Case I with a = 0.005 and P = 0.1 
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The analytic PE transmission loss TLp^ is thus 
" “201ogio|P| 

= lOlogio-^-2Ologio|/3o/ro(0 - aQ\nr)\ - 201og 

Figure 10 displays TLp^ at a depth of 1000m. Both transmission loss results 
show that the azimuth variation of sound speed produces lower loss at large 
azimuth angle. 

The difference between the analytic PE transmission loss TLp^ and the 
numerical NM transmission loss is shown in Figure 11. Since the 

modal transmission loss is a function of amplitude alone this difference has 
the same shape as the relative error in the amplitude of the envelope function. 
The difference in transmission loss is everywhere less than 1 dB. 

B. CASE II 

In the second case, we examined a sound speed field that varies in range 
and azimuth, closely representing an eddy or ring structure in the real ocean. 

Figures 12 and 13 show the amplitude and phase, respectively, of the 
envelope function from the normal mode model. Note that azimuth 

v'ariation of the amplitude is larger than range variation. Figures 14 and 15 
are the amplitude and phase of the envelope function Upg based on the 
analytic solution to the wide-angle PE. Figure 16 shows the relative 
difference of the amplitude. The percent difference is everywhere less than 
2%. The difference between the NM and PE phases is shown in Figure 17. 
The difference of phase is everywhere less than one degree. 



7T 

Sin z • 
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Figure 9. at a Depth of 1000 m for Case I with a = 0.005 and p = 0.1 




Figure 10. TLpg at a Depth of 1000 m for Case I with a = 0.005 and p = 0.1 
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Figure 11. Difference in TL in dB for Case I with a = 0.005 and P = 0.1 
The NM and PE transmission losses, TLj^j^ and TLpg, at a depth of 1000 m 
are shown in Figures 18 and 19, respectively. They are nearly the same as can 
be seen in Figure 20 where the error is everywhere less than 0.2 dB. 
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Figure 12. Amplitude of for Case II with a = 4x10 ^ and j3 = 0.1 
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Figure 14. Amplitude of Up^ for Case II with «= 4x10 ^ and P = 0.1 
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Figure 16. Amplitude Difference in Percent for Case II with a = 4x10 and^ = 

0.1 
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Figure 19. TLpg at a depth of 1000 m for Case II with a = 4x10 ^ and /3 = 0.1 
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Figure 20. Difference in TL in dB for Case II with a = 4x10 ^ and P = 0.1 
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V. DISCUSSION 



In this chapter we address two issues. The first is a quantification of 
horizontal refraction. The second is an assessment of model accuracy. 

A. HORIZONTAL REFRACTION 

Critical to this accuracy test is the use of sound speed fields that cause 
significant horizontal sound refraction. To confirm that the selected sound 
speed fields do cause significant horizontal refraction, we compare the 
corresponding Nx2D NM solutions (the calculation is divided into N vertical 
slices and a 2D NM model is used in each radial direction) with 3D NM 
solutions. 

The envelope functions of the N x 2D solutions have been calculated 
from (A.25) in the appendix. The amplitude of the NM envelope function is 
the same for both the Nx2D and 3D solutions in Case I and Case II. However, 
the phases are different. Figure 21 shows the phase of the Nx2D envelope 
function for Case I. The phase is essentially constant at -45°. In Case II, the 
phase of the envelope function of the Nx2D solution is also constant and its 
value is also -45°. 

Now let us compare the phases calculated by the Nx2D method with the 
phases of the 3D solutions for Case 1 and Case II, as displayed in Figures 4 and 
13, respectively. Both cases show a 13°-15° phase difference over a range of 
100 km, implying that the azimuthal sound speed variation used is large 
enough to induce sound propagation out of the vertical plane. 
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Figure 21. Phase of for Nx2D for Case I with a = 0.005 and P = 0.1 

B. ACCURACY 
1. Phase 

In this section w'e discuss the differences between the analytic PE and 
numerical NM solutions. In Case I, the phase difference between the analytic 
PE and numerical NM solutions is less than 2.8 degree everywhere. The Case 
II phase differences between the numerical NM and analytic PE solutions, as 
shown in Figure 17, is everywhere less than 1 degree. A quantification of the 
phase differences for the two runs is given in Table 2. 



TABLE 2. PHASE DIFFERENCE 



Case 


PE 


Normal Mode NM 


Difference IPE-NMi 




min 


max 


min 


max 


max 


Case I 


-45.15 


-60.35 


-45.15 


-58.92 


2.8 


Case II 


-45.00 


-59.04 


-45.00 


-58.94 


0.43 
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2. Amplitude 

Results are shown in Table 3. Envelope amplitude of NM and PE 
generally agree well with the largest relative difference near 5%. The largest 
difference occurs in a region that has the largest azimuth gradient in sound 
speed. 



TABLE 3. AMPLITUDE DIFFERENCE 



Case 


Amplitude of NM 


Amplitude of PE 


Error Value 
Difference 


Relative 
Difference [%] 




min 


max 


min 


max 


max 


max 


Case I 


0.798 


4.22 


0.817 


4.27 


0.038 


4.8 


Case II 


0.851 


4.28 


0.841 


4.28 


0.076 


1.8 



3. Transmission Loss (TL) 



The difference in transmission loss between the analytic PE and 
numerical NM solutions is shown in Table 4. The difference for both cases is 
everywhere less than half a decibel. 



TABLE 4. TRANSMISSION LOSS DIFFERENCE 



Case 


TL of NM 


TL of PE 


Difference 




min 


max 


min 


max 


max 


Case I 


28.84 


63.31 


28.85 


63.73 


0.43 


Case II 


28.73 


62.75 


28.73 


62.85 


0.15 
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VI. CONCLUSIONS 



In this work we have tested a three-dimensional numerical acoustic 
coupled mode model. We have sought to describe its accuracy in regard to 
the modeling of horizontal refraction. In order to test its accuracy, we 
compared the coupled mode model results with two different analytic 
solutions to the parabolic wave equation. 

We have concentrated on the accuracy of the envelope function and 
transmission loss calculation. For the sound speed fields described in Table 
1, we found that the phase error of the slowly varying envelope function is 
lower than 2 degrees and that the envelope amplitude agrees to within 5% or 
better. Our test cases have indicated that the normal mode model agrees 
closely with the analytic solutions. 

We have only used depth-in\^ariant sound speed fields for this test. As a 
result, the accuracy in modeling mode-mode interactions is not tested here. 
Future tests should be directed at examining the accuracy of calculating 
mode-coupling effects using depth varying sound speed fields. 
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APPENDIX 



A. THREE-DIMENSIONAL COUPLED NORMAL MODE METHODS 
1. Overview 

a. Three-Dimensional (3D) Helmholtz Equation. 

The acoustic pressure P(f, t) in the ocean is governed by the wave 

equation; 






1 d'^P{r,t) _ 

C^(r) dt^ 



(AT) 



where C(r) represents the speed of sound propagation. The cylindrical 

coordinates are r = (r, 0, z) where r is range, d is the azimuthal angle 

2 

measured positive counterclockwise, -z is depth, and V is the Laplacian 
operator in cylindrical coordinates. 



„2 1 a 1 a“ a^ 



(A.2) 



Inclusion of an acoustic source q{r,t), modifies Equation (AT) to 

= (^- 3 ) 

C of 

For an acoustic point source located at Pq with time-harmonic circular 
frequency w, Equation (A.3) becomes 

= -4KS{f - (AA) 

If we let 
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(A.5) 



p(r,l)=0(r)c 

where (P(r) is time-independent, then a Helmholtz equation results: 



V^0(j) + k^{r)0{f) = -47r5{f -r^) 



(A.6) 



where the acoustic wavenumber is k(r) = 
Equation is 



TV 

c{r) ■ 



For y Q, the Helmholtz 



V^0(f) + k^(f)0(y) = 0 . 



(A.6a) 



b. Coupled Normal- Mode Solution 

Adiabatic normal mode theory was applied to a range-dependent 
medium by Pierce [1965] and Milder [1969]. The pressure field was expressed 
as a linear combination of the local modes whose coefficients were obtained 
from a system of decoupled ordinary differential equations. The physics of 
coupled normal mode theory includes non-adiabaticity and results in a 
coupled system of differential equations. The theory used in the 
development of the Chiu-Ehret 3D coupled normal mode model follows. 

The first step in the development is the expansion of the time- 
independent pressure component in terms of the local normal modes, i.e., 

0(') = ^'LPn(r,»)Zn(^-r,e) (A.7) 

where 

is the n*^ local normal modes at point (r,0). 
is the mode amplitude function at point (r,0). 

c. The Local Normal Modes 

The local normal modes obey the following equation: 
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dz^ 



• + 






z„=o 



(A.8) 



th 



where is the horizontal wavenumber associated with the n mode at each 
horizontal location. 

The idealized upper and lower boundary conditions are 

Z„(z = 0;r,6>) = 0. (A.9) 



^Z„(z = -H;r,0) = O. 
oz 



(A.10) 



where H is the ocean depth. The local normal modes can be normalized 
according to the orthonormal condition 



H 



~ ^mn- 



(A.n) 



d. The Mode Amplitude Function 

To obtain the governing equation for P^, substitution of (A. 7) in 
(A. 6a) is necessary. The use of (A.8) and the farfield approximation 



1.2 ii 



mm TH m m 



(A.12) 



following the substitution results in the following equation: 



I 



m 






d^P 1 d^P ^ 
dr^ 3^2 



+ 2 



dZ^ ^ 1 dP^ dZ^ 



dr dr r2 dd d6 



+R 



^ 1 3^Z„ 



m 



dr^ 



r^ dd^ 






(A.13) 



The next step is to multiply (A.13) by Z„ and then integrate over z, i.e.. 
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z 



m 



j 



^3^ 1 3^ 2p ' 

3r2 aa- 



dz + 2}z„l 



dP^, dz„i 1 dP,„ dz 



m '-'^m 



df dr dd dd 



dz 



'j^rn^n 



^d^Z^ 1 d-Z, 



\ 



+ 



m 



dr^ dd^ 



dz 

y j 



0 . 



(A.14) 



A subsequent application of the orthonormality condition (A. 11) on (A.14) 
gives 



2 13 

+ k„ + 



2 A 



^02 

' r2 302 



= -Z 



m 



2JZ„ 



dz. 



m 



dr 



dz 



dr r 






az. 



rn 



dz 

de j 



dR 



m 



dd 



Pmj^n 



d^Z^ . 1 d^Z^ 



dr^ de- 



dz 



Defining coupling coefficients as 



(A. 15) 



y„„ = 2jz„^rfz 

/Sm« = -JZ„ 

r •' 



dZ 



de 



dz 



^mn ~ j 



^d^Z^ 1 d^Z. 



m 



dr^ r^ de- 



dz, 



the coupled system (A.15) can be recast as 



^02 

3r2 



7 2 10 

+ kn + 



2 \ 



r2 302 






m 



' A « 1A^ 

/““31 '^'”%3«J 



4- B P 



(A16) 



e. The Envelope Functions 

Following the work of Chiu and Ehret (1990), the mode 
amplitude function P^ can be separated into a slowly varying envelope 
function and a rapidly varying component i.e.. 
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<Pn = \lK{r, 0 )dr. 

A substitution of (A. 17) into (A. 16) yields 



(A17) 

(A18) 



3 ^ ,2 13 ^'' 

3 r^ r^de^j 



Uj*- = -E 



m 



fy — + B -— 1+B U e’^ 

Imn-^ ^ Hmn ^ ^mn 

V or r dPy 



which can be expanded as 



f3' 1 3^ 1 

1 


•n, ^ 

4_ t 9 ]> 4_ J 


dk„ i d^(pn 1 


2^ 


1 2/ a</»„ a 


dO^ j 


dr 


r^ 80^ V y 




30 30_ 






= -I 



m 



Ymn + Pmn + {^YrnnKn + ^Pmn — ^mn 

dr r dO J \ r d6 

JO. 



Hi 



(A19) 
(A. 20) 



Dividing (A.20) by i2k^e " gives the governing equation for the envelope: 



Dr 



^ d 2 1 



0 r^ r^ dd^ 



\ 


fd J 


L/,+ 

J 


V + En 
^or J 



U„ + F„ l^U„ = XG„„ ■ vu„ + 

^ m m 

(A. 21) 



where 



D„ 



-I 



2 ^:. 



d^„ . _^J_f d(p„ ^ 

\d0 j 



^r dO^ 



k^r 



do 



^mn 



2kr, 



rmnr + PmnOy^'^^ 
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V 



(m n) 

(m = «) 



d ^ 13; 
— r + -—0 
dr r do 



-1 



H 



2k 



ytnn^m 



Pfnn 



m 



mn 



n L 



do 



-iB 



mn 






2k. 



■B 



nn 



and f and 0 are the unit directional vectors in r and 0, respectively. 

/. Transmission Loss TL 

The mode amplitude function is related to the envelope 

function and the phase by (A. 17). The local normal modes are 

solutions to the eigenvalue-eigenfunction problem (A.8)-(A.ll). Once and 

are computed, pressure can be calculated from the product of and Z^. 

2 

At each point (r, 0, z), 0 is 

<X>-(r,3,z) = [Re(0)f -t[lm(<P)f (A.22) 



where 

Re[0{r,e,z)] = -Z[Re{P„)]z„(z) 

m 

and 

Im(0) = X[lm(P^)]-Z^(z) 

m 

In (A. 22), Re and Im are used to denote the real and imaginary 
parts, respectively. Transmission loss can be computed as 

TL = -101ogio<2>Ve,z). (A.23) 
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2. Numerical Solution 

a. Numerical Model 

To obtain a normal mode solution we must solve the equations 
for the eigenvalues and eigenmodes (A.8)-(A.ll) and the equation for the 
modulation envelope (A. 21). These equations can be solved numerically. 

b. The Eigenvalues and the Eigenmodes 

Equation (A. 8) is approximated using central finite differencing. 
This approximation casts (A.8)-(A.10) into a matrix algebra problem. The 
eigenvalues and the eigenmodes of which can be determined using the 
iterativ'e QR method (Acton, 1970). 

c. The Iterative Method 

To obtain the envelope function L/„, we need to solve equation 
(A. 21) which can be rearranged to form a set of first order partial differential 
equations (PDE) with smaller terms put on the right-hand side. We can solve 
for L/ji in an iterative fashion. The iterative equation is 



dr 



+ £, 






H, 



m 



lA 

r ^ 



U' = - F - — U’ - D 



dr 



,2 



1 

r^ dO^ 



U 



/■-I 



IC 



mu ^ u 



(A. 24) 



m 



where U „ is the solution at the iteration for mode n. This set of first order 
differential equations is integrated using a Runge-Kutta method of order 5 
and 6 (Acton, 1970). 

3. The N X 2D Method 



In the N x 2D method, the sound channel is divided into N vertical 
slices and the two-dimensional normal mode solution is solved for each slice. 
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It combines the results in each vertical plane to construct a 3D field. This 
method provides an approximation to a 3D solution for a 3D sound field. 

Thus, in each slice the envelope function U ,, is independent of 
azimuthal variation effects. With the azimuthal terms ignored, the governing 
equation becomes 



dr 



U 



i -1 



m 






f-1 



m 



where 



D„ = 






Ik 



E„ = 



n 

1 dk 



n 



2k ij dr 



r = — !_ y 
^uiii I mil''' 

AKij 

V = — r 
dr 



H 



inn 



2k, 



’iynm^in ‘ ^^inn)^ 



i( 0HI ■ 'Pii) 






■B 



inn 



(A25) 



(m * n) 
{m = n) 



B. THREE-DIMENSIONAL PARABOLIC WAVE EQUATION 
1. Parabolic Wave Equation (PE) 

The following work was developed by Tapper! (1977) and Lee (1988). 
In cylindrical coordinates, the Helmholtz Equation (A. 6a) becomes 



1 a0 1 d^O 
+ + — 

dr^ r dr r^ dO^ 



d^0 .2 2 

+ — + kgn 0{r,6 ,z) = 0 
dz 



(A.26a) 



where 
/Cq = w/Cq 
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Cq is a reference sound speed 

XV is the harrrionic source frequency 

n is the index of refraction 
n = n{r,e, 2 ) = Cq /C{r,9, z). 

The PE approximation assumes 

0{r,e,z) = u{r,e,z)v{r) 



where v{r) is rapidly varying and u{r, 9, 2 ) is a modulation. 
By substituting (A. 26b) into (A. 26a), one obtains 



d^u 

-U 


' 1 ^ 2dv' 


du 1 d^ti d^u 

-j- - -j- 


w 


, r V dr , 


dr r^ d9^ dz^ 



2 2 

+ kQii u 





d^v 


1 dll 


V + 


- — - + — 


dr^ 


r dr 

J 



u 



If xj is determined by 



d^v 1 dv 

:r+ + 

r dr 



kg V = 0 



then II must satisfy 



dr 2 



+ 



1 ^ 2 di' 
r dr , 



dll 



1 d^ii 



dr d9^ 



d^ii 

dz^ 



A'o - 1)1/ = 0. 



The solution to (A. 27) is 



v[r) = nlikQr) ~ 



1 



r > 

2 


2 


r ~ 

. , n 


^nkQr^ 


exp 


• ‘o’- ■ J 

v ^ y - 



where Hq is a Hankel function of zeroth order. 

By substituting (A. 29) into (A. 28), one gets 



(A.26b) 



= 0. (A. 26c) 



(A.27) 



(A.28) 



(A.29) 
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d 1 , 2 / 

~rr ~ + 2 ->^2 + ^0 (” ~ ^) 

dr 



dr dz^ r^ 39 ^ 



u = 0 



(A.30) 



which can be factored into 



d d 

— + ikQ - IRqQ — + ikQ + ikQQ 
,dr A or 



u = 0 



(A.31) 



where 



Q = 



1 + 






kQ dz 



+ 



ky 



02 

dO^ 



2 

2 



Q is called the square root operator. If one considers only outgoing waves, the 
equation to be solved is 



' 3 •, •, 

- + ikQ- ikoQ 

dr J 



u = 0. 



(A.32) 



Let 



so that 



1 

Q = [l + X + y]2 (A.33) 




y = 



_1 ^ 

ky dd^ 



and then expand Q in a Taylor series: 

Q = [i + x + y]^ «i + -x--x2+ly+.... 



8 



(A.34) 
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Inclusion of the first four terms in the expansion leads to a "wide angle" 
version of the PE approximation. Substituting (A.34) into (A. 32) gives 



du _ 
dr 2 




kQ dz"^ 




1 

kQ dz^ 



2 



ky d# 



(A.35) 



2. The Envelope Function of the Parabolic Equation 

Here we present the analytical results of Seigmann et al., (1990). Their 
development uses the pressure release surface and rigid bottom boundary 
conditions, i.e., 

i/(r,0,2 = 0) =0 

|^(r,0,2= H) =0. (A.36) 

dz 



A solution to equation (A.35) for a depth-independent sound speed field and 
subject to the boundary conditions stated above can take the following form: 



u{r,9,z) = L/^(r,0)sin()'y2)exp 







2 






2 


ikQV 


Yi 




1.1 


y) 




2 






4 




-I, 



(A.37) 



where 







n 



H' 



H is the ocean depth, and U ^ is the envelope function of the parabolic 
solution. 

By substituting (A.37) into (A.35), one obtains the equation governing 
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Uf,ir, 0). 



dUf>(r,6) _ 

dr 2 



(2 


r 

, 1 


in - 1) 


1 +- 


V / 


K) 



^21 



7/ 





y) 


4 9 

1 a^ 


U'oJ 


{korf^e^^ 



(A.38) 



Letting 



UAr,6) = 



py,^i 

where both A (r,0) and 0(r,0) are real quantities, (A.38) becomes 



(A. 39) 



dA 



dr 



V . d© 

+ lA.^ = 

P dr 2 



(n^ - 1) 



1 

1 +- 
2 



r 

IL 

L^OJ 



1/2 1 
- hr -1 - - 

A \ ' 0 



y) 


4^ 


U'oJ 


j 



Ik^r 



d Ap a A., a© 

+ 2i _ + lA, 



dQ^ 



de de 



d^Q 

d^ 



lae J 



A, 



(A. 40) 



A separation of the real and imaginary parts of (A.40) yields 

dA, 



q; 1 dApde A q2q 
+ :r + ^ = 0 



dr kQr^ d9 dO Ik^p- d^ 



(A.41) 



[„^ - 1) 



- 


r > 


2" 






4 




lA 

2 






1/2 1 
4^ ^2 


7; 


2 a© 1 

— _j_ 


'a©' 


l^'oj 




^'0 

V. J 


*■'0 (korf 


. de. 



1 1 

[kQrf Ap de^ 

(A.42) 



Given A^ and Q, a pressure field is determined, and the index of refraction can 
be calculated from (A.42). The resultant sound speed field can then be used 
in the normal mode model. 
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